#' ----------------------------------------------------
#' C7 paper on individual legislative efforts
#' Landtag BW - EJPR reviews additional analyses
#' Second R&R
#' Lion Behrens
#' ----------------------------------------------------

#' ------------------------------
# 1. supplementary analyses -----
#' ------------------------------

  #' ---------------------------------------
  # 1.1 do newbies behave differently? -----
  #' ---------------------------------------
  
    #' --------------------------------------------
    # 1.1.1 comparison of amendment behavior ------
    #' --------------------------------------------
  
      par(mfrow=c(1,2))
    
      # n_amends 
      hist(termxmp$n_amends[termxmp$newbie==0], prob = TRUE, breaks=20,
           col = "grey",
           main = "Returners", 
           xlab = "Number of Submitted Amendments")
      par(new = TRUE)
      boxplot(termxmp$n_amends[termxmp$newbie==0], horizontal = TRUE, axes = FALSE,
              col = rgb(0, 0.8, 1, alpha = 0.5))
      
      hist(termxmp$n_amends[termxmp$newbie==1], prob = TRUE, breaks=20,
           col = "grey",
           main = "New MPs", 
           xlab = "Number of Submitted Amendments", 
           ylab = "")
      par(new = TRUE)
      boxplot(termxmp$n_amends[termxmp$newbie==1], horizontal = TRUE, axes = FALSE,
              col = rgb(0, 0.8, 1, alpha = 0.5))
   
    
      # n_articles
      hist(termxmp$n_articles[termxmp$newbie==0], prob = TRUE, breaks=20,
           col = "grey",
           main = "Returners", 
           xlab = "Number of Submitted Article Changes")
      par(new = TRUE)
      boxplot(termxmp$n_articles[termxmp$newbie==0], horizontal = TRUE, axes = FALSE,
              col = rgb(0, 0.8, 1, alpha = 0.5))
      
      hist(termxmp$n_articles[termxmp$newbie==1], prob = TRUE, breaks=20,
           col = "grey",
           main = "New MPs", 
           xlab = "Number of Submitted Article Changes", 
           ylab = "")
      par(new = TRUE)
      boxplot(termxmp$n_articles[termxmp$newbie==1], horizontal = TRUE, axes = FALSE,
              col = rgb(0, 0.8, 1, alpha = 0.5))
    
      
    #' -----------------------------------------------------------------
    # 1.1.2 tracing new MPs' amendments over the legislative term ------
    #' -----------------------------------------------------------------
      
      # timing of amendments 
      billxmp$id <- 1:nrow(billxmp)
      amend_timing <- billxmp %>% 
        group_by(id) %>% 
        summarize(idDrucksache = first(idDrucksache),
                  datedist_days = first(datedist_days),
                  days_passed = first(days_passed), 
                  newbie = newbie, 
                  amend = amend
        ) %>% 
        filter(newbie == 1, 
               amend == 1)
      
      par(mfrow=c(1,1))
      hist(amend_timing$days_passed, breaks=200, xaxt="n", yaxt="n",
           main="Timing of Amendments for New MPs", xlab="Number of Days since Last Election")
      axis(side=1, at=seq(0,1800, 100), labels=seq(0,1800, 100))
      axis(side=2, at=seq(0,14,2), labels=seq(0,14,2))
      
    
      
    
    #' --------------------------------------------
    # 1.1.3 subset analysis for final models ------
    #' --------------------------------------------
      
      # define data and controls
      billxmp_opp <- billxmp[which(billxmp$opposition==1),]
      billxmp_gov <- billxmp[which(billxmp$opposition==0),]
      dat_gov <- billxmp_gov[-which(billxmp_gov$candidatures=="99"),]
      dat_opp <- billxmp_opp[-which(billxmp_opp$candidatures=="99"),]
      
      controls <- ~. + weiblich + age + experience + committee_member + position_highest + Mandat + urban.district + log(district_magnitude) + vote_distance + 
        log(factionsize) + log(datedist) + oppchair_committee + log(n_bills_field) + log(bill_length) + cap1_field + (1|idDrucksache) + (1|mp_name)
      
      # models for newbies
      m1_newbie <- glmer(update(amend ~ candidatures, controls), 
                  family="binomial", nAGQ=0, data = dat_opp[dat_opp$newbie==1,])
      performance::icc(m1_newbie, by_group=T)
      r.squaredGLMM(m1_newbie)
      
      m2_newbie <- glmer(update(amend ~ media, controls), 
                  family="binomial", nAGQ=0, data = dat_opp[dat_opp$newbie==1,])
      performance::icc(m2_newbie, by_group=T)
      r.squaredGLMM(m2_newbie)
      
      m3_newbie <- glmer.nb(update(n_articles ~ candidatures, controls), nAGQ=0, data = dat_opp[dat_opp$newbie==1,])
      performance::icc(m3_newbie, by_group=T)
      r.squaredGLMM(m3_newbie)
      
      m4_newbie <- glmer.nb(update(n_articles ~ media, controls), nAGQ=0, data = dat_opp[dat_opp$newbie==1,])
      performance::icc(m4_newbie, by_group=T)
      r.squaredGLMM(m4_newbie)
      
      # models for non-newbies
      m1_NOnewbie <- glmer(update(amend ~ candidatures, controls), 
                         family="binomial", nAGQ=0, data = dat_opp[dat_opp$newbie==0,])
      performance::icc(m1_NOnewbie, by_group=T)
      r.squaredGLMM(m1_NOnewbie)
      
      m2_NOnewbie <- glmer(update(amend ~ media, controls), 
                         family="binomial", nAGQ=0, data = dat_opp[dat_opp$newbie==0,])
      performance::icc(m2_NOnewbie, by_group=T)
      r.squaredGLMM(m2_NOnewbie)
      
      m3_NOnewbie <- glmer.nb(update(n_articles ~ candidatures, controls), nAGQ=0, data = dat_opp[dat_opp$newbie==0,])
      performance::icc(m3_NOnewbie, by_group=T)
      r.squaredGLMM(m3_NOnewbie)
      
      m4_NOnewbie <- glmer.nb(update(n_articles ~ media, controls), nAGQ=0, data = dat_opp[dat_opp$newbie==0,])
      performance::icc(m4_NOnewbie, by_group=T)
      r.squaredGLMM(m4_NOnewbie)
      
      stargazer(m1_NOnewbie, m2_NOnewbie, m3_NOnewbie, m4_NOnewbie, m1_newbie, m2_newbie, m3_newbie, m4_newbie, 
                covariate.labels = c("Progressive Ambition", "Static Ambition", "Social Media Presence", "Female", "Age",
                                     "Seniority", "Committee Member", "Party Group Leadership", "Parliamentary Leadership", "District Winner", "Closeness of Electoral Race",
                                     "Urban", "Log (District Magnitude)", "Log (Party Group Size)", "Log (Date Distance)", "Opposition Chair Committee", "Log(Number of Bills in Policy Field)", "Log (Bill Length)"), 
                star.cutoffs = c(.05, .01, .001)) 
      
      

    #' --------------------------------------------------------------------
    # 1.2 restrict models of Table 4 to last year of legislative term -----
    #' --------------------------------------------------------------------
    
      # timing of bills 
      legis_timing <- billxmp %>% 
        group_by(idDrucksache) %>% 
        summarize(idDrucksache = first(idDrucksache),
                  datedist_days = first(datedist_days),
                  days_passed = first(days_passed)
        )
      
      par(mfrow=c(1,1))
      hist(legis_timing$days_passed, breaks=200, xaxt="n", yaxt="n",
           main="Timing of Bills", xlab="Number of Days since Last Election")
      axis(side=1, at=seq(0,1800, 100), labels=seq(0,1800, 100))
      axis(side=2, at=seq(0,14,2), labels=seq(0,14,2))
      
      threshold <- min(legis_timing$days_passed[which(legis_timing$datedist_days < 366)])
      rect(threshold, 0, 1800, 12, density = 50, angle = 45,
           col = "lightblue", border=NA)
      text(x = 1625, y = 12.3, "Last year of legislative period", cex = 0.82)
      
      
      # second: re-fit models 1-4 only for activity in last year 
      # variable datedist is left out and cap1_field is left out for convergence
      dat_opp$leadership <- 0
      dat_opp$leadership[dat_opp$position_highest == "faction leadership" | dat_opp$position_highest == "parliamentary leadership"] <- 1
      controls <- ~. + experience + newbie + committee_member + leadership + vote_distance +
        log(factionsize) + oppchair_committee + log(n_bills_field) + log(bill_length) + (1|idDrucksache) + (1|mp_name)
      
      m1_lastyear_simple <- glmer(amend ~ candidatures + (1|idDrucksache) + (1|mp_name), 
                  family="binomial", nAGQ=0, data = dat_opp[dat_opp$datedist_days<366,])
      m1_lastyear_controls <- glmer(update(amend ~ candidatures, controls), 
                  family="binomial", nAGQ=0, data = dat_opp[dat_opp$datedist_days<366,])
      
      m2_lastyear_simple <- glmer(amend ~ media + (1|idDrucksache) + (1|mp_name), 
                                    family="binomial", nAGQ=0, data = dat_opp[dat_opp$datedist_days<366,])
      m2_lastyear_controls <- glmer(update(amend ~ media, controls), 
                  family="binomial", nAGQ=0, data = dat_opp[dat_opp$datedist_days<366,])
      
      m3_lastyear_simple <- glmer.nb(n_articles ~ candidatures + (1|idDrucksache) + (1|mp_name), data = dat_opp[dat_opp$datedist_days<366,])
      m3_lastyear_controls <- glmer.nb(update(n_articles ~ candidatures, controls), nAGQ=0, data = dat_opp[dat_opp$datedist_days<366,])
      
      m4_lastyear_simple <- glmer.nb(n_articles ~ media + (1|idDrucksache) + (1|mp_name), nAGQ=0, data = dat_opp[dat_opp$datedist_days<366,])
      m4_lastyear_controls <- glmer.nb(update(n_articles ~ media, controls), nAGQ=0, data = dat_opp[dat_opp$datedist_days<366,])
      
      stargazer(m1_lastyear_simple, m1_lastyear_controls, m2_lastyear_simple, m2_lastyear_controls,
                m3_lastyear_simple, m3_lastyear_controls, m4_lastyear_simple, m4_lastyear_controls,
                covariate.labels = c("Progressive Ambition", "Static Ambition", "Seniority", "First Term", "Committee Member", "Leadership Position", 
                                     "Closeness of Electoral Race", "Log (Party Group Size)", "Opposition Chair Committee", "Log(Number of Bills in Policy Field)", "Log(Bill Length)", "Social Media Presence"), 
                star.cutoffs = c(.05, .01, .001)) 
      
      
      # vary the threshold for all four models, plot coefficients 
      coefs_prog <- coefs_static <- coefs_media <- 
        as.data.frame(matrix(NA, nrow=length(seq(540, 220, -20))*2, ncol=4))
      colnames(coefs_prog) <- colnames(coefs_static) <- colnames(coefs_media) <- 
        c("model", "n_days", "coef", "std.err")
      coefs_prog$model <- coefs_static$model <- 
        c(rep("m1", nrow(coefs_prog)/2),
          rep("m3", nrow(coefs_prog)/2))
      coefs_media$model <- 
        c(rep("m2", nrow(coefs_prog)/2),
          rep("m4", nrow(coefs_prog)/2))
      coefs_prog$n_days <- coefs_static$n_days <- coefs_media$n_days <-
        rep(seq(540, 220, -20),2)
      
      for (n_days in seq(540, 220, -20)) {
        
        # fit models
        m1 <- glmer(update(amend ~ candidatures, controls), 
                             family="binomial", nAGQ=0, data = dat_opp[dat_opp$datedist_days<n_days,])
        m2 <- glmer(update(amend ~ media, controls), 
                             family="binomial", nAGQ=0, data = dat_opp[dat_opp$datedist_days<n_days,])
        m3 <- glmer.nb(update(n_articles ~ candidatures, controls), nAGQ=0, data = dat_opp[dat_opp$datedist_days<n_days,])
        m4 <- glmer.nb(update(n_articles ~ media, controls), nAGQ=0, data = dat_opp[dat_opp$datedist_days<n_days,])
        
        # store coefficients, m1
        coefs_prog[which(coefs_prog$model=="m1" & coefs_prog$n_days ==n_days), "coef"] <-
          fixef(m1)["candidaturesprog"]
        coefs_prog[which(coefs_prog$model=="m1" & coefs_prog$n_days ==n_days), "std.err"] <-
          sqrt(diag(vcov(m1)))[2]
        
        coefs_static[which(coefs_static$model=="m1" & coefs_static$n_days ==n_days), "coef"] <-
          fixef(m1)["candidaturesstatic"]
        coefs_static[which(coefs_static$model=="m1" & coefs_static$n_days ==n_days), "std.err"] <-
          sqrt(diag(vcov(m1)))[3]
        
        # store coefficients, m2
        coefs_media[which(coefs_media$model=="m2" & coefs_media$n_days ==n_days), "coef"] <-
          fixef(m2)["media"]
        coefs_media[which(coefs_media$model=="m2" & coefs_media$n_days ==n_days), "std.err"] <-
          sqrt(diag(vcov(m2)))[2]
        
        # store coefficients, m3
        coefs_prog[which(coefs_prog$model=="m3" & coefs_prog$n_days ==n_days), "coef"] <-
          fixef(m3)["candidaturesprog"]
        coefs_prog[which(coefs_prog$model=="m3" & coefs_prog$n_days ==n_days), "std.err"] <-
          sqrt(diag(vcov(m3)))[2]
        
        coefs_static[which(coefs_static$model=="m3" & coefs_static$n_days ==n_days), "coef"] <-
          fixef(m3)["candidaturesstatic"]
        coefs_static[which(coefs_static$model=="m3" & coefs_static$n_days ==n_days), "std.err"] <-
          sqrt(diag(vcov(m3)))[3]
        
        # store coefficients, m4
        coefs_media[which(coefs_media$model=="m4" & coefs_media$n_days ==n_days), "coef"] <-
          fixef(m4)["media"]
        coefs_media[which(coefs_media$model=="m4" & coefs_media$n_days ==n_days), "std.err"] <-
          sqrt(diag(vcov(m4)))[2]
        
        print(n_days)
        
      }

    
      
        
        # -------------------------------------------------
        # coefficient for progressive ambition (M1, M3)
        # -------------------------------------------------
        
          # m1
          plot(coefs_prog[coefs_prog$model=="m1", "n_days"], 
               coefs_prog[coefs_prog$model=="m1", "coef"], 
               ylim=c(-4,8), xaxt="n",
               pch=19, col="lightblue", lwd=2,
               xlab="Number of days left in legislative period", 
               ylab="Coefficient estimate",
               main="Coefficient for Progressive Ambition"
          )
          axis(1, at=sort(coefs_prog[coefs_prog$model=="m1", "n_days"])+2.5, 
               labels=sort(coefs_prog[coefs_prog$model=="m1", "n_days"]))
          
          abline(h=0, col="red", lty=2, lwd=2)
          segments(coefs_prog[coefs_prog$model=="m1", "n_days"], 
                   coefs_prog[coefs_prog$model=="m1", "coef"] + 1.96 *coefs_prog[coefs_prog$model=="m1", "std.err"], 
                   coefs_prog[coefs_prog$model=="m1", "n_days"], 
                   coefs_prog[coefs_prog$model=="m1", "coef"] - 1.96 *coefs_prog[coefs_prog$model=="m1", "std.err"]
          )
          
          # m3
          points(coefs_prog[coefs_prog$model=="m3", "n_days"]+5, 
                 coefs_prog[coefs_prog$model=="m3", "coef"], 
                 ylim=c(-4,8), pch=19, col="orange", lwd=2,
                 xlab="Number of days left in legislative period", 
                 ylab="Coefficient estimate"
          )
          
          segments(coefs_prog[coefs_prog$model=="m3", "n_days"]+5, 
                   coefs_prog[coefs_prog$model=="m3", "coef"] + 1.96 *coefs_prog[coefs_prog$model=="m3", "std.err"], 
                   coefs_prog[coefs_prog$model=="m3", "n_days"]+5, 
                   coefs_prog[coefs_prog$model=="m3", "coef"] - 1.96 *coefs_prog[coefs_prog$model=="m3", "std.err"]
          )
          
          legend(490, 7.5, legend=c("Model 1", "Model 3"),
                 col=c("lightblue", "orange"), lwd=2,
                 box.lty=0)
        
      
      
          # -------------------------------------------------
          # coefficient for static ambition (M1, M3)
          # -------------------------------------------------
          
            # m1
            plot(coefs_static[coefs_static$model=="m1", "n_days"], 
                 coefs_static[coefs_static$model=="m1", "coef"], 
                 ylim=c(-2,4), xaxt="n",
                 pch=19, col="lightblue", lwd=2,
                 xlab="Number of days left in legislative period", 
                 ylab="Coefficient estimate",
                 main="Coefficient for Static Ambition"
            )
            axis(1, at=sort(coefs_static[coefs_static$model=="m1", "n_days"])+2.5, 
                 labels=sort(coefs_static[coefs_static$model=="m1", "n_days"]))
            
            abline(h=0, col="red", lty=2, lwd=2)
            segments(coefs_static[coefs_static$model=="m1", "n_days"], 
                     coefs_static[coefs_static$model=="m1", "coef"] + 1.96 *coefs_static[coefs_static$model=="m1", "std.err"], 
                     coefs_static[coefs_static$model=="m1", "n_days"], 
                     coefs_static[coefs_static$model=="m1", "coef"] - 1.96 *coefs_static[coefs_static$model=="m1", "std.err"]
            )
            
            # m3
            points(coefs_static[coefs_static$model=="m3", "n_days"]+5, 
                   coefs_static[coefs_static$model=="m3", "coef"], 
                   ylim=c(-2,4), pch=19, col="orange", lwd=2,
                   xlab="Number of days left in legislative period", 
                   ylab="Coefficient estimate"
            )
            
            segments(coefs_static[coefs_static$model=="m3", "n_days"]+5, 
                     coefs_static[coefs_static$model=="m3", "coef"] + 1.96 *coefs_static[coefs_static$model=="m3", "std.err"], 
                     coefs_static[coefs_static$model=="m3", "n_days"]+5, 
                     coefs_static[coefs_static$model=="m3", "coef"] - 1.96 *coefs_static[coefs_static$model=="m3", "std.err"]
            )
            
            legend(490, 3.5, legend=c("Model 1", "Model 3"),
                   col=c("lightblue", "orange"), lwd=2,
                   box.lty=0)
        
        
            # -------------------------------------------------
            # coefficient for social media presence (M2, M4)
            # -------------------------------------------------
            
              # m2
              plot(coefs_media[coefs_media$model=="m2", "n_days"], 
                   coefs_media[coefs_media$model=="m2", "coef"], 
                   ylim=c(-1,2), xaxt="n",
                   pch=19, col="lightblue", lwd=2,
                   xlab="Number of days left in legislative period", 
                   ylab="Coefficient estimate",
                   main="Coefficient for Social Media Presence"
              )
              axis(1, at=sort(coefs_media[coefs_media$model=="m2", "n_days"])+2.5, 
                   labels=sort(coefs_media[coefs_media$model=="m2", "n_days"]))
              
              abline(h=0, col="red", lty=2, lwd=2)
              segments(coefs_media[coefs_media$model=="m2", "n_days"], 
                       coefs_media[coefs_media$model=="m2", "coef"] + 1.96 *coefs_media[coefs_media$model=="m2", "std.err"], 
                       coefs_media[coefs_media$model=="m2", "n_days"], 
                       coefs_media[coefs_media$model=="m2", "coef"] - 1.96 *coefs_media[coefs_media$model=="m2", "std.err"]
              )
              
              # m4
              points(coefs_media[coefs_media$model=="m4", "n_days"]+5, 
                     coefs_media[coefs_media$model=="m4", "coef"], 
                     ylim=c(-1,2), pch=19, col="orange", lwd=2,
                     xlab="Number of days left in legislative period", 
                     ylab="Coefficient estimate"
              )
              
              segments(coefs_media[coefs_media$model=="m4", "n_days"]+5, 
                       coefs_media[coefs_media$model=="m4", "coef"] + 1.96 *coefs_media[coefs_media$model=="m2", "std.err"], 
                       coefs_media[coefs_media$model=="m4", "n_days"]+5, 
                       coefs_media[coefs_media$model=="m4", "coef"] - 1.96 *coefs_media[coefs_media$model=="m2", "std.err"]
              )
              
              legend(490, 2, legend=c("Model 2", "Model 4"),
                     col=c("lightblue", "orange"), lwd=2,
                     box.lty=0)
        
        
      
    #' -------------------------------------------------------------
    # 1.3 investigate vote motivation (closeness of race) ----------
    #' -------------------------------------------------------------
    
       #' -------------------------------
       # 1.3.1 descriptive analysis -----
       #' -------------------------------
       
        legis_activity <- billxmp %>% 
          group_by(mp_name, term) %>% 
          summarize(vote_distance = first(vote_distance),
                    amend = mean(amend)*100,
                    n_articles = sum(n_articles)
          )
        legis_activity$vote_distance_cat <- NA
        legis_activity$vote_distance_cat[legis_activity$vote_distance < 45] <- 1
        legis_activity$vote_distance_cat[legis_activity$vote_distance < 35] <- 2
        legis_activity$vote_distance_cat[legis_activity$vote_distance < 30] <- 3
        legis_activity$vote_distance_cat[legis_activity$vote_distance < 25] <- 4
        legis_activity$vote_distance_cat[legis_activity$vote_distance < 20] <- 5
        legis_activity$vote_distance_cat[legis_activity$vote_distance < 15] <- 6
        legis_activity$vote_distance_cat[legis_activity$vote_distance < 10] <- 7
        legis_activity$vote_distance_cat[legis_activity$vote_distance < 5] <- 8
        
        par(mfrow=c(1,1))
        # y = amend
        boxplot(legis_activity$amend ~ legis_activity$vote_distance_cat, xaxt="n",
                xlab="Closeness of Electoral Race / Distance to Competitor", ylab="Percentage of Bills Amended", 
                main="Absence of Vote-Seeking Behavior in Amendment Efforts")      
        axis(side=1, at=1:8, labels=c("35-45%", "30-35%", "25-30%", "20-25%", 
                                      "15-20%", "10-15%", "5-10%", "0-5%"))   
        
        # y = n_articles
        boxplot(legis_activity$n_articles ~ legis_activity$vote_distance_cat, xaxt="n",
                xlab="Closeness of Electoral Race / Distance to Competitor", ylab="Number of Submitted Article Changes", 
                main="Absence of Vote-Seeking Behavior in Amendment Efforts")      
        axis(side=1, at=1:8, labels=c("35-45%", "30-35%", "25-30%", "20-25%", 
                                      "15-20%", "10-15%", "5-10%", "0-5%"))   
              
        
       #' ---------------------------------------------------------------
       # 1.3.2 refit M1-M4 with different dummy operationalizations -----
       #' ---------------------------------------------------------------
        
        coefs_vote_distance <- 
          as.data.frame(matrix(NA, nrow=length(seq(2,15,1))*4, ncol=4))
        colnames(coefs_vote_distance) <- c("model", "perc", "coef", "std.err")
        coefs_vote_distance$model <- 
          c(rep("m1", nrow(coefs_vote_distance)/4),
            rep("m2", nrow(coefs_vote_distance)/4), 
            rep("m3", nrow(coefs_vote_distance)/4), 
            rep("m4", nrow(coefs_vote_distance)/4) 
            )
        coefs_vote_distance$perc <- rep(seq(2,15,1),2)
        
        for (perc in seq(2,15,1)) {
          
          # define dummy vote distance variable
          dat_opp$vote_distanceCAT <- dat_opp$vote_distance
          dat_opp$vote_distanceCAT[dat_opp$vote_distance<perc] <- 0
          dat_opp$vote_distanceCAT[dat_opp$vote_distance>perc] <- 1
          
          controls <- ~. + weiblich + age + experience + newbie + committee_member + position_highest + Mandat + urban.district + log(district_magnitude) + vote_distanceCAT + 
            log(factionsize) + log(datedist) + oppchair_committee + log(n_bills_field) + log(bill_length) + cap1_field + (1|idDrucksache) + (1|mp_name)
          
          # re-fit models
          m1 <- glmer(update(amend ~ candidatures, controls), 
                      family="binomial", nAGQ=0, data = dat_opp)
          m2 <- glmer(update(amend ~ media, controls), 
                      family="binomial", nAGQ=0, data = dat_opp)
          m3 <- glmer.nb(update(n_articles ~ candidatures, controls), nAGQ=0, data = dat_opp)
          m4 <- glmer.nb(update(n_articles ~ media, controls), nAGQ=0, data = dat_opp)
          
          # store coefficients
          coefs_vote_distance[which(coefs_vote_distance$model=="m1" & coefs_vote_distance$perc ==perc), "coef"] <-
            fixef(m1)["vote_distanceCAT"]
          coefs_vote_distance[which(coefs_vote_distance$model=="m1" & coefs_vote_distance$perc ==perc), "std.err"] <-
            sqrt(diag(vcov(m1)))[13]
          
          coefs_vote_distance[which(coefs_vote_distance$model=="m2" & coefs_vote_distance$perc ==perc), "coef"] <-
            fixef(m2)["vote_distanceCAT"]
          coefs_vote_distance[which(coefs_vote_distance$model=="m2" & coefs_vote_distance$perc ==perc), "std.err"] <-
            sqrt(diag(vcov(m2)))[12]
          
          coefs_vote_distance[which(coefs_vote_distance$model=="m3" & coefs_vote_distance$perc ==perc), "coef"] <-
            fixef(m3)["vote_distanceCAT"]
          coefs_vote_distance[which(coefs_vote_distance$model=="m3" & coefs_vote_distance$perc ==perc), "std.err"] <-
            sqrt(diag(vcov(m3)))[13]
          
          coefs_vote_distance[which(coefs_vote_distance$model=="m4" & coefs_vote_distance$perc ==perc), "coef"] <-
            fixef(m4)["vote_distanceCAT"]
          coefs_vote_distance[which(coefs_vote_distance$model=="m4" & coefs_vote_distance$perc ==perc), "std.err"] <-
            sqrt(diag(vcov(m4)))[12]
          
          print(perc)
          
        }
          
          
        # -------------------------------------------------
        # coefficient for vote distance
        # -------------------------------------------------
        
        # m1
        plot(coefs_vote_distance[coefs_vote_distance$model=="m1", "perc"], 
             coefs_vote_distance[coefs_vote_distance$model=="m1", "coef"], 
             ylim=c(-2,5), xaxt="n",
             pch=19, col="lightblue", lwd=2,
             xlab="Distance Between Candidates That Is Used As a Threshold in %", 
             ylab="Coefficient estimate",
             main="Coefficient for Closeness of the Electoral Race (Dummy)"
        )
        axis(1, at=sort(coefs_vote_distance[coefs_vote_distance$model=="m1", "perc"]), 
             labels=sort(coefs_vote_distance[coefs_vote_distance$model=="m1", "perc"]))
        
        abline(h=0, col="red", lty=2, lwd=2)
        segments(coefs_vote_distance[coefs_vote_distance$model=="m1", "perc"], 
                 coefs_vote_distance[coefs_vote_distance$model=="m1", "coef"] + 1.96 *coefs_vote_distance[coefs_vote_distance$model=="m1", "std.err"], 
                 coefs_vote_distance[coefs_vote_distance$model=="m1", "perc"], 
                 coefs_vote_distance[coefs_vote_distance$model=="m1", "coef"] - 1.96 *coefs_vote_distance[coefs_vote_distance$model=="m1", "std.err"]
        )
        
        # m2
        points(coefs_vote_distance[coefs_vote_distance$model=="m2", "perc"]+0.1, 
               coefs_vote_distance[coefs_vote_distance$model=="m2", "coef"], 
               ylim=c(-1,2), pch=19, col="orange", lwd=2,
               xlab="Distance Between Candidates that is Used as a Threshold", 
               ylab="Coefficient estimate"
        )
        
        segments(coefs_vote_distance[coefs_vote_distance$model=="m2", "perc"]+0.1, 
                 coefs_vote_distance[coefs_vote_distance$model=="m2", "coef"] + 1.96 *coefs_vote_distance[coefs_vote_distance$model=="m2", "std.err"], 
                 coefs_vote_distance[coefs_vote_distance$model=="m2", "perc"]+0.1, 
                 coefs_vote_distance[coefs_vote_distance$model=="m2", "coef"] - 1.96 *coefs_vote_distance[coefs_vote_distance$model=="m2", "std.err"]
        )
        
        
        # m3
        points(coefs_vote_distance[coefs_vote_distance$model=="m3", "perc"]+0.2, 
               coefs_vote_distance[coefs_vote_distance$model=="m3", "coef"], 
               ylim=c(-1,2), pch=19, col="darkgreen", lwd=2,
               xlab="Distance Between Candidates That Is Used As a Threshold in %", 
               ylab="Coefficient estimate"
        )
        
        segments(coefs_vote_distance[coefs_vote_distance$model=="m3", "perc"]+0.3, 
                 coefs_vote_distance[coefs_vote_distance$model=="m3", "coef"] + 1.96 *coefs_vote_distance[coefs_vote_distance$model=="m3", "std.err"], 
                 coefs_vote_distance[coefs_vote_distance$model=="m3", "perc"]+0.3, 
                 coefs_vote_distance[coefs_vote_distance$model=="m3", "coef"] - 1.96 *coefs_vote_distance[coefs_vote_distance$model=="m3", "std.err"]
        )
        
        # m4
        points(coefs_vote_distance[coefs_vote_distance$model=="m4", "perc"]+0.3, 
               coefs_vote_distance[coefs_vote_distance$model=="m4", "coef"], 
               ylim=c(-1,2), pch=19, col="purple", lwd=2,
               xlab="Distance Between Candidates that is Used as a Threshold in %", 
               ylab="Coefficient estimate"
        )
        
        segments(coefs_vote_distance[coefs_vote_distance$model=="m4", "perc"]+0.2, 
                 coefs_vote_distance[coefs_vote_distance$model=="m4", "coef"] + 1.96 *coefs_vote_distance[coefs_vote_distance$model=="m4", "std.err"], 
                 coefs_vote_distance[coefs_vote_distance$model=="m4", "perc"]+0.2, 
                 coefs_vote_distance[coefs_vote_distance$model=="m4", "coef"] - 1.96 *coefs_vote_distance[coefs_vote_distance$model=="m4", "std.err"]
        )
        
        
        
        legend(13, 4, legend=c("Model 1", "Model 2", "Model 3", "Model 4"),
               col=c("lightblue", "orange", "darkgreen", "purple"), lwd=2,
               box.lty=0)
        
        
      
       
        
              
  
         
      
      
      
      
      



